The computation of generalized percolation critical 
: polynomials by the deletion-contraction algorithm 
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Abstract. Although every exactly known bond percolation critical threshold is the 
root in [0, 1] of a lattice-dependent polynomial, it has recently been shown that the 
notion of a critical polynomial can be extended to any periodic lattice. The polynomial 
is computed on a finite subgraph, called the base, of an infinite lattice. For any problem 
with exactly known solution, the prediction of the bond threshold is always correct, 
regardless of the base chosen. For unsolved problems, the polynomial is referred to as 
\ the generalized critical polynomial and provides an approximation that becomes more 

accurate with increasing number of bonds in the base, appearing to approach the exact 
answer. The polynomials are computed using the deletion-contraction algorithm, which 
quickly becomes intractable by hand for more than about 18 bonds. Here, I present 
■ generalized critical polynomials calculated with a computer program for bases of up 

£h ', to 36 bonds for all the Archimedean lattices, except the kagome which was considered 

in an earlier work. The polynomial estimates are generally within 10 -5 to 10~ 7 of the 
numerical values, but the prediction for the (4, 8 2 ) lattice, though not exact, is not 
ruled out by simulations. 
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1. Introduction 

Percolation [T] is one of the simplest random processes taking place on a lattice. 
Nonetheless, since its introduction over fifty-five years ago [2], it has continued to provide 
physicists and mathematicians with an array of fascinating problems (see [21 El [5l EJ [7_1 [8] 
for a small sampling) and has inspired many new mathematical techniques [91 [10], EH]- 
Given an infinite lattice, L, we declare each edge to be open with probability p and 
closed with probability 1 — p. When p is small, L will be sparsely populated by small 
clusters of open bonds. When p is near 1, we will have an infinite open mass with small 
pockets of closed bonds. In between this regime lies the critical threshold, p c , which 
marks the transition from the unconnected phase to the phase containing an infinite open 
cluster. One of the most challenging problems in the field is the analytic determination 
of critical thresholds. Outside of one dimension and a narrow class of two-dimensional 
lattices [T2"l |13"1 H3] , exact results remain elusive. In some cases, mathematically rigorous 
confidence intervals [15] and bounds [Ml [17] have been proved and which are continually 
improving, but, for most lattices, critical probabilities are known only numerically 
[T8l [T9l l20l I2T1 [22] . All solved two-dimensional lattices are formed from self-dual 3- 
uniform hypergraphs, such as the one shown in Figure UK, where the shaded triangle in 
Figure [lb permits any configuration of sites and bonds as long as they lie within the 
three boundary vertices. For such graphs, the critical condition is given by [131 E3J E] 

P(A,B,C) = P(A,B,C) (1) 

where P(A, B, C) is the probability that all three corners are connected and P(A, B, C) 
is the probability that none are connected. For bond percolation, for example, this 
condition allows one to find the inhomogeneous critical surface by assigning each bond 
in the shaded triangle a different probability. The homogeneous bond threshold is given 
by a polynomial in the probability p, of degree equal to the number of bonds in the 
triangle. In previous work [211 1213 [26] it was shown that one can define a critical 
polynomial on any lattice that agrees with ([T]) for self-dual 3-uniform lattices, and 
provides an accurate approximation for unsolved problems. This polynomial is a kind 
of graph invariant [27] with similarity to the Tutte polynomial |28] in that it may be 
computed by the deletion-contraction algorithm. In fact, deletion-contraction is how 
the generalized critical polynomial is defined, as we will see in the following section. An 
infinite lattice may be partitioned by a finite subgraph which is tiled in a regular way 
to give the full lattice. The critical polynomial is defined on this subgraph, called the 
base, and is thus a property of a finite graph and its embedding into the infinite lattice. 
Using the example of the kagome lattice, it was demonstrated in [27] that polynomials 
computed on bases of increasing numbers of bonds provide better approximations to 
the critical thresholds, with bases of 36 bonds making predictions within 10 -7 of the 
numerically determined critical threshold. 

In this paper, I use the computer program used in [27] for the kagome lattice to 
find polynomials for the remaining unsolved Archimedean lattices (Fig. [2]), providing 
further evidence for the conjecture that the estimates made by these critical thresholds 
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Figure 1. a) A self-dual 3- uniform hypergraph. The shaded triangles can be any 
configuration of sites and bonds lying within the three boundary vertices; b) the critical 
point is given by ([!]). 



approach the exact answer as the size of the base increases. I begin by reviewing the 
definition of the polynomial in the next section. In section [31 I describe in detail the 
operation of the program to automate its computation. The subsequent sections are 
devoted to reporting polynomials for the Archimedean lattices partitioned into different 
bases. 

2. Deletion-contraction algorithm 

Consider the hexagonal lattice of Figure [2b, with assignment of probabilities on the unit 
cell shown in Figure [3b,. The condition ([T|) then yields the critical surface 

H (p, r, s) = prs — pr — ps — rs + 1 = . (2) 

As with all surfaces determined by (pQ), this expression is first-order in all its arguments. 
Consider the (4, 8 2 ) lattice of Figure [2k, with the assignment of probabilities in Figure 
[3b. This lattice has an unsolved bond threshold as it does not fall into the class that may 
be found with ([[]). However, we generalize the notion of the critical polynomial on this 
lattice by assuming that its critical surface has the first-order property. Next, we note 
that contracting its p-bond by setting p = 1 yields the martini-A lattice [12]. This lattice 
is formed by substituting the generator in Figure [3b for the shaded triangles in Figure 
[lb, and thus its critical surface is known exactly. It has previously been reported in 
many places [251 EHl ED] so here we just denote it A(r, s, t, u, v). Deleting the p-bond by 
setting p = gives the hexagonal lattice with some bonds doubled in series. Reduction 
to these two cases and the imposition of the first-order property leaves us with only 
one choice for the critical surface of the (4, 8 2 ) lattice, given by the deletion-contraction 
formula, 

FE(p, r, s, t, u, v) = pA{r, s, t, u, v) + (1 — p)H(s, ur, tv) . (3) 
Expanding this, and setting all probabilities equal gives the homogeneous polynomial, 
1 - 4p 3 - 2p A + 6p 5 - 2p 6 = , (4) 
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(a) (b) (c) 

Figure 3. Assignments of probabilities on a) the hexagonal lattice, b) the (4,8 2 ) 
lattice, c) the martini-A lattice. 

with solution on [0, 1], p c = 0.67683519... . The simulation result reported by Parviainen 
|20j is / um = 0.67680232(63) (the brackets indicate the standard error on the last digits), 
which, although the difference is only 3.3 x 10~ 5 , easily rules out our estimate. However, 
we may extend the inhomogeneous probabilities over two unit cells. That is, we may 
consider a larger base. This is shown in Figure [T2k where each of the twelve bonds 
should be understood to have a different probability, and the shapes on the external 
vertices indicate how the base is embedded in the infinite lattice. Now when deleting 
and contracting a bond, we will not immediately see solvable lattices. Nevertheless, we 
can recursively apply the deletion-contraction algorithm until known lattices appear. 
The result of this is a binary tree of graphs, the size of which increases exponentially 
with the number of bonds in the base. For this reason, a computer is needed to handle 
large bases, but even then the maximum size achieved in this work is only 36 bonds. 
However, this is generally sufficient to provide estimates to thresholds within 10 -7 of 
the numerical values. 

The 12 th -order polynomial for the (4,8) base in Figure [T2"a is [26] . 

1 - 4/ - 16/ + 12/ + 22/ + 16/ - 70/° + 48/ 1 - 10/ 2 = 0, (5) 

with solution p c = 0.67678736..., differing from the numerical value by 1.5 x 10~ 5 
and cutting the error in half from the 6-bond base. At first glance, it appears that 
the generalized critical polynomial depends on the order in which bonds are chosen 
in the deletion-contraction algorithm. However, in [27j, an argument was given that 
the generalized critical surface, and therefore the polynomial, is independent of the 
bond order. In the computer implementation, the answers are tested by performing the 
calculations several times using random bond orders. 

3. Implementation 

The computer implementation of this algorithm is written in C++. As input, it takes a 
collection of bonds and a set of vertices, along with the number, N e , of external vertices. 
Each bond and vertex is assigned a number, and the problem is fully specified once the 
two end vertices of each bond are given and identifications are made between external 
vertices to indicate how the base is embedded into the infinite lattice. Two vertices that 
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are identified are called "equivalent" . For a given lattice L, the solve routine works as 
follows: 

(i) Choose a bond, b, on which to perform the deletion and contraction. A minimal set 
of requirements is used here. Of course, the same bond is chosen for both deletion 
and contraction. 

(ii) Delete b to form lattice L , and contract the bond to form L\. 

(iii) Simplify L and L\ by removing or combining any superfluous bonds that appeared 
as a result of deletion or contraction. 

(iv) Check if L$ is one of the known lattices. If so, assign it the proper name and stop. 
If not, recursively call the solve routine for this lattice. 

(v) Check if L\ is one of the known lattices. If so, assign it the proper name and stop. 
If not, recursively call the solve routine for this lattice. 

This process produces a binary tree of lattices, with branches that terminate on identified 
graphs. The run-time is related to the number of lattices that appear in this tree. 
However, the population is determined by the speed with which known lattices appear, 
and only loosely related to the number of bonds, n, in the base. As different lattices and 
bases can have very different connectivity properties, the run-time may vary considerably 
for problems with the same n. But even for a given lattice, the choice of bond, b, at 
each step can significantly change the speed with which the solution is found, although 
the final result is independent of these choices. However, adding a bond to a base 
essentially doubles the number of lattices that will appear in the tree, and thus the 
number of operations increases as roughly 2 n for a given problem. 

3.1. Setup 

Figure H] shows an example of a lattice partitioned into a base with nine bonds and 
nine vertices, eight of which are external. The embedding of the base is specified by 
identifying vertices that are equivalent in the infinite lattice and these identifications are 
illustrated by matching shapes in Figure H^i. The numbering of the bonds and internal 
vertices is completely arbitrary. However, the external vertices are always assigned the 
first N e numbers and their ordering is important. Note that we only consider connected 
bases that can thus be contained inside a simple boundary loop that runs through all 
the external vertices in the obvious way (the dotted curve in Figure UK) . The vertices 
are numbered counter-clockwise around this loop. When we delete bonds and simplify 
lattices, this ordering will be crucial in recognizing equivalent external gaps. 

3.2. Deletion 

When a bond is deleted, it is simply removed from the list. Some simplification of 
the resulting lattice may be necessary. For example, we may create bonds doubled in 
series, like bonds p\ and p 4 , and p 2 and p 3 , of Figure [5}3. These are replaced by single 
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(a) 



(b) 



Figure 4. a) A nine-bond base. The shapes indicate external vertex identifications, 
b) the lattice resulting from tiling this base. 



Pi 





(b) 

l-(l-p 4Pi )(l-p 3 P 2 ) 
(d) 



Figure 5. Deletion of bond in a) puts bonds 1 and 4, and 2 and 3, in series, as 
depicted in b); c) simplifying these bonds leads to two bonds in parallel; d) the result 
of all the simplifications yields one bond with an effective probability. 



bonds with probabilities p\p± and P2P3, as shown in Figure [5b. Although series bonds 
are usually the only direct result of deletion, Figure [5b indicates how simplifying these 
bonds can indirectly lead to bonds doubled in parallel. These bonds also need to be 
replaced with effective bonds with the appropriate probabilities. For example, in Figure 
[6b, which is the result of contracting bond in Figure [6^, bond p\ is deleted, and p 2 is 
assigned the new probability 

1 - (1 -pi)(l -P2) =Pi +P2 -P1P2 (6) 

which is simply the statement that in order for the effective single bond to be open, pi 
and P2 cannot both be closed. 
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(a) 



(b) 



P,+P 2 -P,P 2 



(P+P 2 -P,P 2 )(P+P 4 -P 3 P 4 ) 



(C) 



(d) 




3.3. Contraction 

When a bond is contracted, it is removed from the base, and its end vertices are merged 
into a single vertex. The direct result of this may be to create bonds doubled in parallel, 
as shown in Figure [6b, but Figure [6fc illustrates how simplifying these parallel bonds 
can indirectly lead to bonds doubled in series. Further complications are possible 
here. Consider the base and lattice in Figure HI Contracting the p§ bond in Figure 
UK necessitates the merger of vertices 6 and 7. However, examination of the embedding 
in Figure Hb reveals that this gap is equivalent to the gap between vertices 2 and 3, 
and we must merge these vertices as well. The identification of equivalent gaps will be 
described in section [3751 

3.4- Bond selection 

Despite the earlier comments about the impact of bond selection on run time, I have not 
made any serious effort to optimize the bond choice for efficiency. Nevertheless, there 
is a variety of reasons a bond may be rejected in the present scheme: 

(i) Deleting it will disconnect the base. The algorithm is only used on bases with a 
single connected component. Note that, in some cases, contracting a bond may 
also result in a disconnected base. Consider the situation in Figure HJ Contracting 
bond po necessitates the removal of bond p^ or p$, since they will now span the 
same gap and need to be combined into a single bond. However, removing either 
of them disconnects the base. Before settling on a candidate bond, we must first 
examine the simplified lattices that result upon both deletion and contraction. It 
is possible to disconnect the lattice without disconnecting the base, with the result 
being a lattice strip, as shown in Figure [71 

(ii) It is a "supporting" bond, that is, its endpoints are equivalent external vertices 
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P, 



P» P 2 P 4 



(a) 



(b) 



(c) 



Figure 7. a) a six-bond base for the square lattice; b) the embedding of this base; c) 
the lattice strip resulting from setting p 4 = 0. 



(e.g., any odd bond in Figure Ek). Contracting such a bond collapses the lattice 
and can result in very complicated situations. It is certainly not impossible to deal 
with these, but I have avoided them for simplicity. The exception to this is the 
case of a lattice strip (Figure [Tfc). As this is essentially a one-dimensional problem, 
it can only be critical when at least one of the supporting bonds has probability 1. 
This means that if there is a supporting bond, pi, the critical surface is of the form 
(1 —pi)f = 0, where / is some function of the other probabilities. Thus, the result 
of contracting this bond gives zero, bringing about a quick identification. 

(iii) It connects two inequivalent external vertices that are not adjacent on the perimeter 
of the base, like bond P2 in Figure [HJ Contracting such a bond does not collapse 
the lattice, but nevertheless leads to complications that I chose to avoid. 

After settling on a set of rules, one may be concerned about the existence of lattices for 
which there is no legal bond choice, and in fact there are such lattices for these rules. 
An example is shown in Figure |9j which appropriately has the appearance of an array 
of stop signs. Other examples are of a similar nature. Although removing bonds 1 to 5 
would clearly disconnect the base, it is not so obvious why p or p 6 are bad choices. The 
reason is that contracting, say, po, makes p<i and p§ into parallel bonds through their 
external connections. But combining these into a single bond necessitates removing one 
of them, which disconnects the base. The bond p$ has the same problem. There are 
several ways around this. One would be to re-partition the lattice into a different base 
that allows a good bond choice; there are clearly better bases to select for this kagome- 
like lattice. Another is to drop the restriction against disconnected bases. However, 
for bases of 36 bonds, the maximum considered here, the appearance of such lattices is 
infrequent enough that I simply discard calculations in which they appear. In order to 
check that the final answer is correct, I run the algorithm several times for a given base 
anyway, as previously mentioned, each with a different labelling of the bonds to ensure 
a different path is taken. 
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Figure 8. a) a base for the (3 3 ,4 2 ) lattice; b) its embedding. The bond P2 is rejected 
by the algorithm because it joins external vertices that are not neighbours on the 
boundary loop of the base. 




(a) (b) 

Figure 9. a) a base for the kagome-like lattice in b). The bond selection algorithm 
laid out in section 13.41 has no legal choices here. 

3.5. Simplification 

During the simplification phase, the algorithm checks the base for extraneous bonds 
or sites, iterating until no further changes need to be made. It seeks bonds doubled 
in series, bonds doubled in parallel, stranded vertices and dead ends. It also checks if 
a bond deletion has resulted in a lattice strip, which changes the bond selection rules 
slightly. 

3.5.1. Series bonds If an internal vertex is incident with only two bonds, then those 
bonds are doubled in series, as in Figure [5b, and they are replaced by a single bond, as 
shown in Figure [5b. It is necessary that we consider only internal vertices here, as an 
external vertex between only two bonds is usually not indicative of bonds doubled in 
series. This can be seen in Figure Hk, where, among other examples, vertex 6 is between 
only bonds p$ and p^, but in the full lattice these bonds are not in series. In some cases, 
an external vertex connected with only one bond may indicate a bond doubled in series, 
which would happen for the lattice in Figure [10] if p were deleted, but I do not replace 
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Figure 10. The "arrow" lattice. Deletion of bond puts bonds 1 and 4 in series, but 
these are not simplified by the algorithm. 

these. 

3.5.2. Parallel bonds If two bonds have the same two end vertices, then they are 
doubled in parallel and one must be removed. This is done by deleting the bond with 
the lower number, and assigning the effective probability to the other. 

It is straightforward to find such parallel cases involving internal vertices, but when 
parallel bonds are between equivalent external vertices, special care is needed. Consider 
Figure ITTa . which shows a three-bond base with external vertex connections that lead 
to the lattice shown in Figure fTTb. Clearly, the bonds po and P2 span equivalent gaps, as 
the resulting lattice is just the square lattice with these two bonds doubled in parallel. 
However, Figure [Tib shows a base with the same configuration of bonds but different 
identifications of external vertices. Tiling this base gives the hexagonal lattice, shown in 
Figure fTTbl in the "brick wall" representation, and now po and P2 do not span equivalent 
gaps, even though each end vertex of p is equivalent with an endvertex of p 2 . It is 
necessary to have some way to discriminate between these cases and it is here that the 
correct numbering of the external vertices becomes important. External vertices that 
neighbour on the boundary loop need not have a bond between them. However, we 
can think of their g ctp clS db directed arc of the loop, (vi,v 2 ), oriented in the clockwise 
direction, with V\ < v 2 . The gap between the last external vertex, v Ne _i, and v is 
written (vn b -i, vo). In Figure fTTb. the bond po spans the gap (2, 3). Vertex 2 is equivalent 
to and vertex 3 is equivalent to 1, and substituting these in the gap (2, 3) gives (0,1), 
which is the correctly oriented gap spanned by bond p . However, in the base of Figure 
UK, the bond p$ bridges the gap (6,7) with vertex 6 equivalent with 3 and 7 with 2. 
Substituting these into (6,7) gives (3,2), an arc directed counterclockwise, as V\ > v 2 , 
which indicates that the gaps (6,7) and (3,7) are equivalent, as can be seen in Figure 
Hb. More plainly, two external gaps are equivalent if the low vertex of one is equivalent 
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Figure 11. The base in a) embeds to give the lattice b), and c) leads to d). In a), the 
bonds po and P2 span equivalent gaps, whereas in c) they do not. 



with the high vertex of the other and vice versa. Note that in Figure Hk, all vertices are 
equivalent so this condition is again met. 

3.5.3. Stranded vertices A stranded vertex is one which is not an endpoint for any 
bond. This can happen for external vertices, as discussed in section 13.21 Stranded 
vertices are easily found and removed. However, if another external vertex, v, was 
equivalent only with the stranded vertex, w, then v is no longer external and must be 
demoted to internal as soon as w is removed. For example, deleting po and p\ in Figure 
[TOb leaves vertex 3 stranded. But vertex 1 was equivalent only with vertex 3 and thus 
the demotion of 1 to internal must accompany the removal of 3. This in turn leaves 1 
a dead-end vertex. 

3.5.4- Dead ends A dead end is an internal vertex that is incident with only one bond. 
This may occur upon demotion of an external vertex, as just described. A dead-end 
bond has no impact on the percolation process and is deleted along with the vertex. 

3.6. Identification 

The algorithm can identify a small number of simple graphs. Although the more lattices 
the algorithm is able to identify, the sooner the process terminates, it can be difficult to 
correctly identify graphs in which not all bonds are equivalent because a bond-matching 
procedure would need to be devised. As such, I only check for the square, triangular, 
hexagonal, and one-dimensional lattices. A lattice with three external vertices, three 
bonds, and one internal vertex is the hexagonal lattice; three external vertices, no 
internal, and three bonds is the triangular lattice; three external vertices and two bonds 
is the square lattice; and a single bond is one-dimensional. Furthermore, as previously 
described, if a strip lattice is collapsed, the result is zero. The capability to identify 
these graphs is sufficient to ensure that the algorithm always succeeds, provided it does 
not run into a lattice of the type in Figure [9j 
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4. Results 



Now I list the results obtained by applying this method to various bases on the 
Archimedean lattices, except the kagome, which was already discussed in [27]. For 
numerical results, we rely heavily on the paper of Parviainen [20], which, with the 
exception of (3, 12 2 ) for which Ding et al. [21] provide a more accurate estimate, is still 
the standard reference for these lattices. 



4.1. (4,8 2 ) lattice 

The single unit cell base result for this lattice was reported in [21], and the 12-bond base 
consisting of two unit cells was studied in [26]. In that work, however, only the base in 
Figure [T2k was used. Another option is to wire the external vertices as in Figure fl~2b . 
but the result of this is the same as that of the 6-bond base. Here, we will label different 
bases as (N x M) where N and M denote the number of unit cells in each direction. 
For example, the (2 x 2) base in Figure [T3h results in the polynomial 

1 - V - 8p 6 + 8p 7 + 16p 9 - 52p w + 16p u - 52p 12 + 296p 13 
+ 160p 14 - 432p 15 - 1142p 16 + 712p 17 + 2436p 18 

- 80p 19 - 7714p 20 + 10520p 21 - 6332p 22 + 1872p 23 - 222p 24 , (7) 

predicting p c = 0.67678965..., which is slightly closer to the numerical value, 
0.67680232(63), than the 12-bond base, but is not a great improvement. Another base 
that might also be called (2 x 2) is shown in Figure [T3~b . which is identical to Figure 
[l~3k except for the identifications of external vertices. This change leads to a different 
polynomial, 

1 - 16p 6 + 4/ - 84p w + 14 4p u + 88p 12 + 176p 13 

- 300p 14 - 144p 15 - 1166p 16 + 1248p 17 + 968p 18 + 4440p 19 

- 16394p 20 + 19392p 21 - 11264j9 22 + 3296p 23 - 390p 24 , (8) 

with p c = 0.67681105..., placing us within 8.73113 x 10~ 6 of the numerical value. 

Moving to bases of 36 bonds, two examples are given in Figures IT4"a and [Tib, though 
many more exist. The first results in the polynomial, 

1 - 6p 4 - 12/ + 12/ + 9p 8 + 24p 9 + 30p 10 - 24p u - 16p 12 

- 108p 13 - 168p 14 + 12p 15 - 612j9 16 + 1428p 17 + 1532p 18 

+ 3144p 19 - 5349p 20 - 15528p 21 - 3246p 22 + 41688p 23 + 40750p 24 



- 77076p 25 - 127878jo 26 + 143284p 27 + 243069p 28 - 234972p 29 

- 607984p 30 + 1414836p 31 - 1365693/ 2 + 758376p 33 - 253242p 34 

+ 47616p 35 - 3898p 36 (9) 

giving p c = 0.6767896635..., again closer to the numerical solution. However, the (3 x 2) 
base in Figure [T4b gives the polynomial 

1 - 12p 7 - 3p 8 - 34p 9 - 12p 10 + 6p n - 6p 12 - 36p 13 
+ 372p 14 + 386p 15 + 378j9 16 - 156p 17 - 1248p 18 - 1800p 19 
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base 



Pc 



(2 x 2)a 


0.67678965... 


(2 x 2)6 


0.67681105... 


(2 x 3) 


0.67678966... 


(3 x 2) 


0.67680215... 



Table 1. Polynomial predictions for various bases of the (4,8 2 ) lattice. 



£>-A 



D— / u 



0; 



Figure 12. Two different embeddings of the (2 x 1) base for the (4, 8 2 ) lattice. 





(a) 



Figure 13. Two different embeddings of the (2 x 2) base for the (4, 8 2 ) lattice. 



5655j9 20 - 
183384j9 25 



4388p 21 + 12528p 22 + 39930p 23 + 62136p 24 



215964p 26 + 345672p z ' + 484599p ze + 48816p 



,27 



2s 



■29 



- 3588878p 30 + 6949548p 31 - 6640005p 32 + 3755348j9 33 

- 1284996p 34 + 247584p 35 - 20 7 28p 36 (10) 

with solution p c = 0.67680215..., a prediction that falls within the standard error of the 
numerical result, and is thus not ruled out. However, it should be stressed that no finite 
base will ever give the exact solution, as demonstrated in 



4.2. (3 3 ,4 2 ) lattice 

Although there is no finite base that will give the exact answer for any of these lattices, 
it is not clear in what manner the base must become infinite. In particular, it may be 
possible that a base need only be infinite in one direction, i.e. a strip. The results of 
this section will rule out this possibility for the (3 3 ,4 2 ) lattice (Figure |2]F). We will start 



Computation of percolation critical polynomials 



15 




(a) 



(b) 



Figure 14. 36-bond bases for the (4,8 2 ) lattice; a) (2 x 3), b) (3 x 2) 



with the progression shown in Figure HU in which the bases consist of unit cells stacked 
in the vertical direction. One unit cell (Figure [T5k) gives the polynomial reported in 

A 



1 - 2p - 2p 2 + 3p 3 - p- 

with solution p c = 0.419308168..., fairly different from Parviainen's 
0.41964191(43). Extending the base to the (2 x 1) configuration (Figure [T5"b) gives 
a polynomial that can be written in the factored form, 

(1 - 2p - 2p 2 + 3p 3 - p 4 )(l - 2p + 2p 2 + p 3 - p A ) = (12) 

and we recognize the first term in brackets as the polynomial (II ip . The second term 
contains no root in [0, 1] and therefore the prediction is the same as for the 5— bond 
base. Similarly, stacking three cells as in Figure [T5b gives the factored form, 

(1 - 2p - 2p 2 + 3p 3 - p A ) (1 - p + p 2 ) (1 - 3p + 2p 2 + Ap 3 - 2p 4 - 2p 5 + p 6 ) (13) 

and once again we have the same prediction. It is probably safe to conjecture that this 
trend continues. The story is different if we extend the base in the horizontal direction. 
The polynomial for the (1 x 2) (Figure [T6b) case was reported in [25], with the result 
p c = 0.419614759... . For (1 x 3) (Figure QSb), we find, 

1 - 3p 2 - 8p 3 - 15p 4 - 3p 5 + 220p 6 + 84p 7 - 2052p 8 
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^num 



,12 



,13 



+ 4698p 9 - 5343j9 10 + 3471j9 n - 1231p i: " + 159p l 
+ 30p 14 - 9j9 15 = (14) 

with solution 0.419650951... . The sequence for the remaining bases is shown in Table 
[2J Clearly, these values are converging to a limit, which appears to be similar to the 
numerical value but not actually correct. In particular, the estimates for n = 5, 6, and 
7 indicate that the first seven digits are 0.4196551, but Parviainen's result rules this 
out easily, as it seems extremely likely that the first five digits are in fact 0.41964. This 
indicates that the (1 x oo) base does not provide the exact answer. We can also rule out 
the possibility that we only need a wider strip, say (2 x oo). If this were exact, then so 
would be the case in which we deleted every horizontal and diagonal bond in the second 
row. Contracting the remaining vertical bonds in that row, we would necessarily find 
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N 


Pc 


1 


0.419308168... 


o 
z 


fl 41 Qfi1 47^Q 


3 


0.419650951... 


4 


0.419654761... 


5 


0.419655145... 


6 


0.419655183... 


7 


0.419655187... 



Table 2. Roots of generalized critical polynomials for the (1 x N) bases of the (3 3 , 4 2 ) 
lattice. These appear to converge to a number ruled out by simulations. 



the same formula for (1 x oo) that we already showed was incorrect. Therefore, for this 
lattice, the only way to get the correct threshold is to consider the (N x M) base in 
which both iV and M go to infinity. As such, we turn to the (2x2), (3x2), and (2 x 3) 
cases. For (2x2), the polynomial is, 

1 - 4/ - 8/ - 4/ + 40/ - 108/ + 372/ + 326/ 

- 2640/ - 4132/° + 40124/ 1 - 101829/ 2 + 145944/ 3 

- 134736/ 4 + 82372/ 5 - 32199/ 6 + 6904/ 7 - 236/ 8 

- 224/ 9 + 36/° = (15) 

with solution p c = 0.4196154184..., barely distinguishable from, but still slightly better 
than, the (1x2) answer. The (3 x 2) base yields the polynomial 

1 - 6p 2 - Up 3 + 6/ + 132/ - 58/ - 204/ - 1275/ 
+ 1272/ + 8514/° + 1836/ 1 - 82380/ 2 + 55062/ 3 
+ 504828/ 4 - 903004/ 5 - 3039471/ 6 + 18382050/ 7 

- 48094255/ 8 + 83086308/ 9 - 104909466/° + 100560770/ 1 

- 74085918/ 2 + 41775876/ 3 - 17645241/ 4 + 5303532/ 5 

- 991518/ 6 + 54936/ 7 + 23244/ 8 - 6048/ 9 + 488/° (16) 

predicting p c = 0.4196154196..., which is hardly different from (2x2). The best estimate 
is given by the (2 x 3) base 

1 - 4p 3 - 21/ - 72/ - 32/ + 42/ + 1017/ + 6160/ 

- 5163/° - 49410/ 1 - 78554/ 2 + 466284/ 3 + 1142880/ 4 

- 6161272/ 5 - 1697280/ 6 + 66144060/ 7 - 221367606/ 8 

+ 428244768pl9 - 574288983/° + 567883664/ 1 - 423773169/ 2 

+ 239027052/ 3 - 100024955/ 4 + 29504220/ 5 - 5317833/ 6 

+ 236622/ 7 + 138600/ 8 - 33642/ 9 + 2625/° (17) 

giving p c = 0.4196453185..., differing from the numerical result by 3.6 x 10~ 6 . 
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(a) 



(b) 



(c) 



Figure 15. (N x 1) bases for the (3 3 ,4 2 ) lattice; a) (1 x 1); b) (2 x 1); c) (3 x 1). 





(a) (b) 



Figure 17. a) (2 x 3) and b) (3 x 2) bases for the (3 3 ,4 2 ) lattice. 
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base p c 



(2 x 2) 


0.4196154184... 


(3 x 2) 


0.4196154196... 


(2 x 3) 


0.4196453185... 



Table 3. Threshold predictions for various (3 3 ,4 2 ) bases. 



4.3. (3,12 2 ) lattice 

The single-cell polynomial for the (3, 12 2 ) lattice has been reported in various places 
[25], but I note it here, 

1 - 3p 4 - 6p 5 + 3p 6 + 15/ - 15p 8 + 4p 9 = 0. (18) 

Oddly enough, this polynomial can be written in factored form, 

(1 + p - 2p 3 + p 4 )(l - p + p 2 + p 3 - 7p 4 + 4p 5 ) = 0. (19) 

The second term in brackets is the one that has the root in [0, 1], p c = 0.74042331..., 
and therefore, for some reason, the prediction for this nine-bond base is given by a fifth 
order polynomial. The numerical value given by Parviainen is p£ um = 0.74042195(80), 
whereas Ding et al. [21] give the more recent result 0.74042077(2). Turning to the 
(2 x 2) case (Figure [T8|) . we have the polynomial, 

1 - 6p s - 24p 9 - 12p 10 + 48p n + 18j9 12 - 120p 13 - 120p 14 
+ 252p 15 + 639p 16 - 192p 17 - 720p 18 + 804p 19 - 1302p 20 

- 1024p 21 + 468p 22 - 14652p 23 + 28 6 35p 24 + 79392p 25 

- 181344p 26 - 88368p 27 + 425721p 28 - 17280p 29 - 996084p 30 
+ 1558008p 31 - 1244067p 32 + 603008p 33 - 180096p 34 

+ 30720p 35 - 2304p 36 = (20) 

which does not factor, and predicts p c = 0.74042099..., a difference of 9.6 x 10~ 7 from 
Parviainen's value, but within 2.2 x 10~ 7 of Ding et al. This also compares favourably 
with the approximation given by Ziff and Gu, p c ~ 0.74042081, based on a numerically 
fit universal condition. For the moment, this is as far as we can go with this lattice. 



44. (4,6,12) lattice 

The single unit cell polynomials for the (4,6,12) and (3 2 ,4,3,4) lattices were given in 
[26] . Because of the size of their bases, we can only extend these estimates to two unit 
cells. Two possible bases for the (4, 6, 12) lattice are shown in Figures [T9k and fT9b . 
These turn out to have the same polynomial, 

1 - 12p 6 - 4p 8 + 16p 9 - Up 10 + 32p u - 106p 12 
+ 72p 13 - 248p 14 + 680p 15 - 220p 16 + 2080p 17 - 2414p 18 
- 1944p 19 - 5573p 20 + 6904jo 21 - 1924jo 22 + 24028p 23 
+ 21869p 24 - 48140p 25 - 201044p 26 + 261640p 27 + 198914j9 28 
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Figure 18. (2 x 2) base for the (3, 12 2 ) lattice. 




(a) (b) 

Figure 19. Two 36-bond bases for the (4,6, 12) lattice. These have the same critical 
polynomial. 



+ 148952p 29 - 2204780p 3U + 4153808p 3i - 3985112p 32 

+ 2288132p 33 - 800476p 34 + 158616p 35 - 13734p 36 (21) 

with solution, p c = 0.69375829..., an incremental improvement over the single-cell value, 
p c = 0.69377849... reported in |26j, but within 2.4 x 10~ 5 Parviainen's numerical result, 
p num = 0.69373383(72) . 

4.5. (3 2 ,4,3,4) lattice 

Using the base in Figure I2"0"a . we find 

1 - 4p 2 - 8p 3 + 20p 5 - 84p 6 + 248/ + 970p 8 
- 3748p 9 - 4042p 10 + 41664p n - 101482p 12 
+ 139104j9 13 - 122372j9 14 + 70604p 15 - 25368p 16 

+ 4484j9 17 + 234j9 18 - 256p 19 + 34p 20 = (22) 

giving p c = 0.41412438... , a difference of 1.3 x 10" 5 from Parviainen's 0.41413743, 
but his standard error is 4.6 x 10~ 7 . The single unit-cell prediction found in [22] was 
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(a) (b) 
Figure 20. 20-bond bases for the (3 2 ,4,3,4) lattice. 



0.414120304..., so we have improved on this estimate somewhat. We can do better with 
the base shown in Figure [20b . for which we find the polynomial 

1 - 8p 3 - 36p 4 - 44p 5 + 88/ + 912p 7 + 266p 8 - 7604p 9 
- 1390p 10 + 71480p n - 194138J9 12 + 277116p 13 - 247748p 14 
+ 142984p 15 - 50264p 16 + 7940p 17 + 1010p 18 - 648p 19 + 82p 20 , (23) 

and p c = 0.41414477..., putting us within 7.3 x 10~ 6 of the numerical result. 



4-6. (3 4 ,6) lattice 

The single-cell polynomial for the (3 4 , 6) lattice was reported in [26] with the 
approximation p c = 0.43437077... . Extending to 30 bonds, as shown in Figures [2Tb 
and b, we have, for both these embeddings, 

1 - 8p 3 - 24p 4 + 2p 5 + 58p 6 - 186p 7 - 2Up 8 + 4100p 9 
+ 7053j9 10 - 41066p n - 90261p 12 + 353308p 13 + 81765lp 14 
- 4135392j9 15 - 883466p 16 + 41247262p 17 - 138495809p 18 
+ 271182796p 19 - 372758548p 20 + 383898312p 21 - 304802391p 22 
+ 188649862p 23 - 90989211p 24 + 33847422p 25 - 9501925p 26 
+ 1938036/ 7 - 268902p 28 + 22392p 29 - 823j9 30 (24) 

with p c = 0.43435240... compared with the numerical result p° um = 0.43430621(50). 
The difference is 4.6 x 10 ~ 5 . 



4.7. (3,4,6,4) lattice 

The single-cell 12-bond polynomial for the (3,4,6,4) lattice was reported in [24, 25], 
which made the prediction p c = 0.524821.... The two 24-bond bases in Figures [22k and 
b yield the polynomial 

1 - 4p 3 - 8p 4 - 8p 5 + 18/ - 20p 7 - 126p 8 + 516p 9 
+ 754j9 10 - 1180J9 11 - 6071p 12 + 976p 13 + 34222p 14 
+ 768j9 15 - 257049p 16 + 652088p 17 - 866150p 18 + 730808p 19 
- 412380p 20 + 155436/ 1 - 37456p 22 + 5168p 23 - 304p 24 (25) 
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(a) (b) 
Figure 21. 30-bond bases for the (3 4 ,6) lattice. These yield the same polynomial. 

0^5 — O 





(a) 



(b) 



Figure 22. Bases for the (3,4,6,4) lattice. These have the same generalized critical 
polynomial. 



lattice 




„num 
Pc 


br ly -Pc um | 


(4,8 2 ) 


0.67680215... 


0.67680232(63) 




(3 3 ,4 2 ) 


0.41964531... 


0.41964191(43) 


3.6 x 10~ 6 


(3, 12 2 ) 


0.74042099... 


0.74042077(2) 


2.2 x 10~ 7 


(4,6,12) 


0.69375829... 


0.69373383(72) 


2.4 x 10~ 5 


(3 2 ,4,3,4) 


0.41414477... 


0.41413743(46) 


7.3 x 10~ 6 


(3 4 ,6) 


0.43435240... 


0.43430621(50) 


4.6 x lO" 5 


(3,4,6,4) 


0.52483166... 


0.52483258(53) 


9.1 x 10~ 7 



Table 4. The best polynomial estimates for the Archimedean lattices along with the 
numerical values. All numerics are from except (3, 12 2 ), which is from [2Tj . 



with root p c = 0.52483166..., differing from Parviainen's 0.52483258(53) by only 
9.1 x 10~ 7 . 
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5. Conclusions 

I have presented the computation of generalized critical polynomials using a computer 
program to perform the deletion-contraction algorithm. The results clearly support the 
conjecture that these polynomials predict thresholds converging to the exact answer 
in the limit of an infinite base. Moreover, this seems to suggest that critical bond 
thresholds for most unsolved problems are not algebraic numbers, since no finite base 
can give the correct threshold unless the single-cell answer is exact. A summary of 
the best results for each lattice is in Table HI Their accuracy ranges from 4.6 x 10 -5 
for the (3 4 , 6) lattice, to the result for the (4, 8 2 ) lattice which is not ruled out by the 
numerical answer. Many paths for future work are available. One might hope that far 
greater accuracy is possible if a method can be found to compute the polynomials more 
efficiently than by using deletion-contraction and indeed this is a subject of current 
study. It should be noted that the generalized critical polynomials presented here are 
related to F.Y. Wu's homogeneity approximation [31] for the q— state Potts model |32j, 
which he has so far been able to apply to get a 6-bond kagome |33j polynomial and, in 
recent work with W. Guo [34J, a 6-bond (4, 8 2 ) polynomial. In particular, his method 
predicts (jlj) in the q — > 1 limit. Recently, the deletion-contraction method was also 
generalized to the q— state Potts model [35J, which allows accurate determination of 
critical points for large bases and arbitrary q, including predictions in the unphysical 
anti-ferromagnetic region. Finally, polynomials can be defined in three dimensions with 
deletion-contraction in exactly the same way as in two, but it remains to be seen if they 
can give any information about the critical point. 

This work was performed under the auspices of the U.S. Department of Energy by 
Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. 
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